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Abstract 

Three-dimensional unsteady RANS computations are 
performed for the flow past a pure heaving/plunging 
rectangular wing to study the effect of reduced frequency 
and heaving amplitude on the thrust generation and 
propulsive efficiency. The numerical solution has been 
obtained by using an implicit Reynolds-averaged Navier- 
Stokes solver IMPRANS that employs finite volume nodal 
point spatial discretization scheme with dual time stepping. 
The results are obtained in the form of aerodynamic 
coefficients, thrust coefficient and propulsion efficiency and 
compared with the available data in the literature. 
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Introduction 

It is well known that pure "plunging" or "heaving" 
wing can produce both lift and thrust with certain 
combinations of plunging amplitude and frequency. 
The origin of thrust generation for an oscillating 
aerofoil was found by Knoller and later independently 
by Betz. The ability of a sinusoidally plunging aerofoil 
or wing can generate thrust. This effect is known as 
the "Knoller-Betz or Katzmayer effect". The Knoller- 
Betz effect was demonstrated in a wind tunnel 
experiment by Katzmayer. Also some experimental 
and computational investigation of the Knoller Betz 
effect was done by Jones et ah Following this, a 
number of studies for theoretical and numerical 
models of oscillating aerofoils were developed by 
Garrick and Lighthill for thin aerofoils oscillating in 
inviscid flow. They significantly overestimated the 
propulsion efficiency (Isogai et al., Ramamurti and 
Sandberg) observed in the low-Reynolds-number 
separated flows in natural flight and in the flight of 
Micro Air Vehicles (Jones et ah). 

The Navier-Stokes simulations were performed by 



Tuncer and Platzer who have predicted more accurate 
flow patterns and propulsive efficiency. Neef and 
Hummel have computed the thrust coefficient values 
for aerofoils and three-dimensional high aspect ratio 
wings by using Euler solver. The reduced frequency is 
limited to about 0.1 to simulate the flow over the 
cruising large birds. In that work, they investigated 
the asymmetric thrust distribution during each period, 
and numerically demonstrated that the additional lift 
has no influence on thrust generation. Jones et ah have 
developed a flapping wing model for testing in a low 
speed wind-tunnel, and compared direct thrust 
measurements with the panel method results for 
several configurations. 

Okamoto and Yasuda experimentally investigated 
aerodynamic characteristics of unsteady flow over the 
wings at low Reynolds number. The aerodynamic 
forces and moment acting on wings of aspect ratio 6 
with heaving and feathering oscillations in a wind 
tunnel were measured at a low Reynolds number less 
than 104. 

Jones et ah have investigated both experimentally and 
numerically for a sinusoidal flapping finite aspect- 
ratio wing with pure heave motion. They used a flat 
plate theory, two and three-dimensional panel codes, 
Euler and Navier-Stokes solvers for their numerical 
analysis. The reduced frequency, mean angle of attack, 
aspect ratio and Reynolds number are varied in their 
computations. 

The present work describes the unsteady flow 
simulation over a rectangular wing moving with pure 
sinusoidal plunge motion for different combinations of 
plunge amplitude and reduced frequency in each case. 

IMPRANS Solver 

An implicit Reynolds-averaged Navier-Stokes solver, 
IMPRANS developed in-house for steady and 
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unsteady compressible viscous flows is employed for 
the present simulations. It solves unsteady RANS 
equations in three-dimensions in a moving domain. 
Dual time stepping approach is used with an implicit 
finite volume nodal point spatial discretization scheme. 
Inviscid flux vectors are calculated by using the flow 
variables at the six neighbouring points of hexahedral 
volume. 

The Reynolds-averaged Navier-Stokes equations for 
three-dimensional unsteady compressible flow in a 
moving domain in non-dimensional conservative form 
are given by 

dU dE dF dG n 

— + — + — + — = 0. (i) 

dt dx dy dz 1 > 

Here, U is the vector of conserved variables, E, F and G 
flux vectors, (x, y, z) is the Cartesian coordinate system 
and t is the time variable. 

An implicit finite volume nodal point scheme with 
dual time stepping approach is employed to solve the 
above governing equations (1). The dual time stepping 
consists of an implicit discretization in real time and 
the marching of solution in a pseudo time to steady 
state at each physical time step. Use of an implicit 
second order accurate backward difference formula 
for discretization in real time and Euler's implicit time 
differencing formula for pseudo time results in the 
following equation 
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Here U m = U (f) = U (m At*) is the solution vector at 
pseudo time level m, A U m = U m+1 - U m is the change in 
U m over the time step At* and At denotes the real or 
physical time step that is required to resolve the 
physical unsteadiness of the flow. The barred 
quantities denote the solution vectors at the previous 
real time levels n and n - 1 whereas R represents the 
spatial operators which give rise to the flux residual 
after a discretization in space. 

This basic equation (2) of the implicit dual time 
stepping technique can be solved at each real time step 
by employing a finite volume nodal point spatial 
discretization. In this approach, the flow variables are 
associated with each mesh point (z, j, k) of the grid and 
the centroids of the eight neighbouring hexahedron 
cells surrounding the nodal point are joined to form 
the control volume Qijk. To facilitate the finite volume 
formulation, the equations are written in integral form 
and the surface integrals are evaluated by summing 
up the contributions due to the flux terms over the six 
faces of the computational cell. Applying integral 



conservative equations to each control volume, 
linearising the changes in flux vectors using Taylor's 
series expansions in time, assuming locally constant 
transport properties, and dropping the superscript m 
we obtain 



A- 6 ^- \AU 
dx 



(3) 
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Here Qijk is the control volume surrounding the nodal 
point (z, k) of the curvilinear grid; A = dEi / dU, B = dFi 
/dU,C = dGi 1 5LI, Er = dEvi/ dlL, Fs = dFvi/ dU y and Gr = 
5Gv3/ dllz are the Jacobian matrices; Ei, Fi and Gi are the 
inviscid flux vectors and Ev, Fv and Gv are the viscous 
flux vectors; Smx, Sm y and Smz are the x, y and z 
components of the surface vector corresponding to the 
m th surface of the control volume. 

It is important to note that the terms containing 
inviscid flux vectors can be calculated by using the 
flow variables at the six neighbouring points and 
Taylor's series expansions can be utilised to discretize 
the derivatives in the viscous flux terms directly in the 
physical plane. The resulting block tridiagonal system 
of equations is solved by using a suitable block 
tridiagonal solution algorithm and proper initial and 
boundary conditions. In order to ensure convergence 
and to suppress oscillations near shock waves, a blend 
of second and fourth order artificial dissipation terms 
is added explicitly. Implicit second order dissipation 
terms are also added to improve the practical stability 
bound of the implicit scheme. The algebraic eddy 
viscosity model due to Baldwin and Lomax is used for 
turbulence closure. 

For a moving body, the equations are solved in the 
inertial frame of reference by employing a grid which 
remains fixed to the body and moves along with it. At 
each real time step t + At, starting from the solution at 
the previous time step t, the solution is marched in 
pseudo time t* using local time stepping. Since the 
choice of physical time step At is no longer limited by 
stability considerations, a much larger time step, with 
a fixed but small number of inner iterations in pseudo 
time, can be used to reduce the undesirably large 
computational time for unsteady flow calculations. 
Based on this dual time stepping method an implicit 
Reynolds averaged Navier-Stokes solver IMPRANS 
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has been developed at CSIR-National Aerospace 
Laboratories for computing a wide variety of two- 
dimensional and three-dimensional unsteady viscous 
compressible flows. This RANS solver has been 
extensively validated for computing unsteady flow 
past pitching aerofoils and wings, plunging aerofoils, 
helicopter rotor blades, wind turbines etc. Here, the 
solver is used for three-dimensional unsteady 
compressible viscous flows over a pure-plunging 
rectangular wing. 

Results 

For all the present simulations, a structured single 
block C-H grid around a NACA 0014 rectangular 
wing of aspect ratio 8 is generated using the 
commercial software package Gridgen V 15.5. The far- 
field boundary is kept at 20 chords away from the 
wing surface. The volume grid used for the present 
computations is shown in Fig. 1(a). While the surface 
grid on the wing is shown in Fig. 1(b). The number of 
grid points is 247 x 65 x 75 in the chord wise, normal 
and span-wise directions respectively. The points are 
clustered properly near the leading and trailing edges 
and near the tip of the wing, the first grid spacing 
normal to the wall being 2.O10 5 c. 




FIG. 1(a) C-H TOPOLOGY GRID AROUND THE WING 




FIG. 1(b) SURFACE GRID ON THE WING 
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The three-dimensional unsteady flow over a pure- 
plunging rectangular wing has been simulated for 
different input parameters using implicit Reynolds 
averaged Navier-Stokes solver IMPRANS. The cases 
considered closely approximate the experimental 
model of Jones et al. Three different plunge amplitudes, 
y = 0.1c, 0.2c and 0.4c with different reduced 
frequencies of k = 0.4, 0.6, 0.8 and 1.0 are considered. 

For all unsteady simulations, steady solution is first 
obtained. After steady state convergence is reached, 
the unsteady flow computation is started with a 
prescribed sinusoidal plunging motion. Five 
consecutive cycles are computed to obtain periodic 
solution. The final cycle results are used to calculate 
time averaged thrust coefficient and propulsive 
efficiency. 

The motion of a single finite aspect ratio wing with a 
pure plunge motion in normal direction is defined by 

y(t*)=haCos(kf) (4) 
and the plunging velocity is given by the following 
expression 

y(t*)=-haksin (kf) (5) 
where non-dimensional heave amplitude, ha= y I c, 
reduced frequency, k = coclll™ and non-dimesional time 
is t* (i.e., pseudo time). 

The mean thrust coefficient and propulsion efficiency 
are computed using the following expressions. 

The mean or time-averaged thrust coefficient is defined 
as 

C t=-C d+(Cd) steady (6) 

where C d is the mean drag coefficient averaged for 
one period of heaving oscillation. (Q) steady IS the steady 
drag of the non-moving wing at its present angle of 
attack. 

The % propulsion efficiency can be calculated from the 
ratio between power output and power input, in this 
case, it is given by 

% Propulsion efficiency, (rjprop) = (C t)/(C v ) (7) 
where mean power input coefficient, C v is calculated 
from the product of lift coefficient, Ci and plunging 
velocity, y(f). 



TABLE 1 AERODYNAMIC PARAMETERS USED FOR HEAVING 
RECTANGULAR WING 



Input flow parameters 




Free stream Mach Number (M-) 


0.3 


Free stream Reynolds Number (Re~ ) 


l.OxlO 6 


Angle of attack (a) 


2° 


Reduced frequency (k) 


0.4, 0.6, 0.8 and 1.0 


Plunge or heave amplitude (y ) 


0.1c, 0.2c and 0.4c 
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The input flow parameters used for the present 
computations are listed in Table 1 for heaving 
rectangular wing. 

Case A (ha = OA, k = OA, 0.6, 0.8 and 1.0) 

In this case, the reduced frequency is varying from 0.4 
to 1.0 in steps of 0.2 with the plunging amplitude of 
0.4. Table 2 shows the time-averaged thrust coefficient 
and propulsion efficiency for the rectangular plunging 
wing. The propulsion efficiency values for k = 0.4 and 
k=1.0 are compared with the available results of Jones 
et al. At the reduced frequency, k = 0.4 the present 
computed value is 67.38% and the Jones et al. values 
are 64.15% (NS (FLOWer)), 67.29% (EULER (FLOWer)) 
and 68.66% (PANEL). Similarly at k = 1.0 the present 
value is 49.73% and the Jones et al. values are 51.09% 
(NS (FLOWer)), 56.01% (EULER (FLOWer)) and 
56.55% (PANEL). The present computed values are in 
good agreement with the results of Jones et al. 

From the Table 2, it can be concluded that the higher 
thrust coefficient of 0.12864 is obtained at higher 
reduced frequency, while on the contrary for higher 
efficiency of 67.38% obtained at lower reduced 
frequency. 

TABLE 2 COMPARISON OF MEAN-THRUST COEFFICIENT AND PROPULSION 
EFFICIENCY 



k 


C t 
(Present) 


°/o 77 prop 

(Present ) 


°/o 7]prop 

(Jones et al. [11]) 


0.4 


0.03545 


67.38 


64.15 (NS (FLOWer)) 
67.29(EULER(FLOWer)) 
68.66 (PANEL) 


0.6 


0.06208 


59.85 




0.8 


0.09379 


54.65 




1.0 


0.12864 


49.73 


51.09 (NS (FLOWer)) 
56.01 (EULER (FLOWer)) 
56.55 (PANEL) 



The variation of lift coefficient with heave distance at 
different reduced frequencies is shown in Fig. 2. 
The lift coefficient values are higher during down- 
stroke than those during up-stroke. The computed 
loops of the aerodynamic coefficients clearly 
demonstrate the hysteretic property existing between 
the up-stroke and down-stroke in all the three cases. 
The variation of thrust coefficient with the heave 
distance at different reduced frequencies is plotted in 
Fig. 3. The thrust coefficient values are smaller during 
the first half of the down-stroke compared to the 
second half of up-stroke and become higher during 
the second half of down-stroke than those during the 
first half of up-stroke. 
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FIG. 2 THE VARIATION OF LIFT COEFFICIENT WITH HEAVE 
DISTANCE FOR THE HEAVING WING AT ha = 0.4 
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FIG. 3 THE VARIATION OF THRUST COEFFICIENT WITH 
HEAVE DISTANCE FOR THE HEAVING WING AT ha = 0.4 




FIG. 4 SURFACE PRESSURE CONTOUR ON LOWER SURFACE 
OF HEAVING WING FOR ONE COMPLETE CYCLE OF 
OSCILLATION AT ha = 0.4 AND k = 0.6 

Case B (ha = 0.1, k = OA, 0.6, 0.8 and 1.0) 

In this case, the reduced frequency is varying from 0.4 
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to 1.0 in steps of 0.2 with heaving amplitude of 0.1. 
Table 3 shows the time-averaged thrust coefficient and 
propulsion efficiency on the rectangular plunging 
wing for different reduced frequencies. From the 
results, it is observed that the higher thrust coefficient 
of 0.008398 is obtained at higher reduced frequency, 
while on the contrary for higher efficiency of 66.81% 
obtained at lower reduced frequency. 

TABLE 3 MEAN THRUST COEFFICIENT AND PROPULSION EFFICIENCY 



3C 



k 


C t 


% f^prop 


0.4 


0.002187 


66.81 


0.6 


0.003926 


60.84 


0.8 


0.005989 


57.31 


1.0 


0.008398 


54.87 



The variation of lift coefficient with heave distance at 
different reduced frequencies is shown in Fig. 5. The 
lift coefficient values are higher during down-stroke 
than those during up-stroke. The variation of thrust 
coefficient with the heave distance at different reduced 
frequencies is plotted in Fig. 6. The thrust coefficient 
values are smaller during the first half of down-stroke 
compared to the second half of up-stroke and become 
higher during the second half of down-stroke than 
those during the first half of up-stroke. 
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FIG. 5 THE VARIATION OF LIFT COEFFICIENT WITH HEAVE 
DISTANCE FOR THE HEAVING WING AT ha = 0.1 
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FIG. 7 SURFACE PRESSURE CONTOUR ON LOWER SURFACE 
OF HEAVING WING FOR ONE COMPLETE CYCLE OF 
OSCILLATION AT ha = 0.1 AND k = 1.0 

Case C (ha = 0.2, k = OA, 0.6, 0.8 and 1.0) 

In this case, the reduced frequency is varying from 0.4 
to 1.0 in steps of 0.2 with the plunging amplitude of 
0.2. Table 4 shows the time-averaged thrust coefficient 
and propulsion efficiency on the plunging rectangular 
wing. The higher thrust coefficient of 0.03338 is 
obtained at higher reduced frequency, while on the 
contrary for the higher efficiency of 67.57% obtained at 
lower reduced frequency. 

TABLE 4 MEAN THRUST COEFFICIENT AND PROPULSION EFFICIENCY 



k 


C t 


% f^prop 


0.4 


0.00888 


67.57 


0.6 


0.01575 


60.80 


0.8 


0.02374 


56.42 


1.0 


0.03338 


53.70 





FIG. 6 THE VARIATION OF THRUST COEFFICIENT WITH 
HEAVE DISTANCE FOR THE HEAVING WING AT ha = 0.1 



FIG. 8 THE VARIATION OF LIFT COEFFICIENT WITH HEAVE 
DISTANCE FOR THE HEAVING WING AT ha = 0.2 
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FIG. 9 THE VARIATION OF THRUST COEFFICIENT WITH 
HEAVE DISTANCE FOR THE HEAVING WING AT ha = 0.2 




pressure plots, it is observed that near the tip of the 
wing initially more variations have been observed at 
the beginning of the cycle (i.e., downward stroke) than 
those at the middle and again more variations at the 
end of the cycle (i.e., upward stroke). 

The comparison of mean thrust coefficient and 
propulsion efficiency at different heave amplitude 'ha 
values are shown in Figs. 11 and 12 respectively, for 
the heaving wing. As 'k' or 'ha increases, the time- 
averaged thrust coefficient increases and propulsion 
efficiency decreases for the range considered. 



ytr - -0.1-114 J. 




FIG. 10 SURFACE PRESSURE CONTOUR ON LOWER SURFACE 
OF HEAVING WING FOR ONE COMPLETE CYCLE OF 
OSCILLATION AT ha = 0.2 AND k = 0.8 

The unsteady variation of lift coefficient with heave 
distance at different reduced frequencies is shown in 
Fig. 8. The lift coefficient values are higher during 
down-stroke than those during up-stroke. The results 
plotted in Fig. 9 show the unsteady variations of thrust 
coefficient with the heave distance at different reduced 
frequencies. The thrust coefficient values are smaller 
during the first half of down-stroke compared to the 
second half of up-stroke and become higher during 
the second half of down-stroke than those during the 
first half of up-stroke. 

For all the three cases, the surface pressure fields are 
plotted in Figs. 4, 7 and 10 on the lower surface of the 
heaving rectangular wing for one complete cycle with 
non-dimensional heave amplitudes of 0.4, 0.1 and 0.2 
at k = 0.6, 1.0 and 0.8 respectively. From the all surface 



FIG. 11 THE COMPARISON OF MEAN THRUST COEFFICIENT 
AT DIFFERENT NON-DIMENSIONAL HEAVE AMPLITUDE 'ha 
VALUES FOR THE HEAVING WING 



FIG. 12 THE COMPARISON OF PROPULSION EFFICEINCY AT 
DIFFERENT NON-DIMENSIONAL HEAVE AMPLITUDE 'ha' 
VALUES FOR THE HEAVING WING 

Conclusions 

The unsteady flow over a rectangular heaving wing 
has been computed by using implicit Reynolds 
averaged Navier-Stokes solver IMPRANS. The effect 
of heaving amplitude and reduced frequency on the 
propulsive efficiency and time-averaged thrust 
coefficient has been studied. The computed results 
agree well with the available data in the literature. 
From the results, it is observed that as reduced 
frequency increases, the time-averaged thrust 
coefficient increases and propulsion efficiency 
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decreases for the range of values considered in the 
present computations. In all the cases, higher thrust 
occurred at higher reduced frequency, while higher 
propulsive efficiency occurred at lower reduced 
frequency. The higher thrust coefficient of 0.12864 is 
obtained at higher reduced frequency of 1.0, while on 
the contrary for higher efficiency of 67.38% obtained at 
lower reduced frequency of 0.4 at the heaving 
amplitude of ha = 0.4 and the higher thrust coefficient 
of 0.008398 is obtained at higher reduced frequency of 
1.0, while for higher efficiency of 66.81% at lower 
reduced frequency of 0.4 the heaving amplitude of 
ha = 0.1. And the higher thrust coefficient of 0.03338 is 
obtained at higher reduced frequency of 1.0, while 
higher efficiency of 67.57% at lower reduced frequency 
of 0.4 the heaving amplitude of ha = 0.2. As heave 
amplitude, ha increases, the time-averaged thrust 
coefficient increases and propulsion efficiency 
decreases for the range considered in the present work. 
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